# Imports
from coopr.pyomo import *
from coopr.opt import SolverFactory
import ReferenceModelBase
import numpy

"""
# Model
model = AbstractModel()

# Sets
model.I = Set() #node
model.J = Set() #node
model.S = Set() #source node
model.D = Set() #demand node

# Data
model.Arc =  Param(model.I, model.J) #arc available
model.Rev =  Param(model.I, model.J) #arc revenue
model.Cost = Param(model.I, model.J) #arc cost
model.B = Param()

# Random
model.ArcDemand = Param(model.I, model.J) #arc demand

# Variables
model.X = Var(model.S, bounds=(0.0, None))
model.Y = Var(model.I, model.J, bounds=(0.0, model.B))
model.Z = Var(model.I, model.J, bounds=(0.0, None))

model.FirstStageProfit = Var()
model.SecondStageProfit = Var()

# Constraints
def vehicle_num_cap_rule(model):
    return sum(model.X[s] for s in model.S) == model.B
model.VehicleNumCapRule = Constraint(rule=vehicle_num_cap_rule)

def source_balance_rule(model, s):
    return sum(model.Y[s,j] for j in model.J if model.Arc[s,j]>=1) == model.X[s]
model.SourceBalanceRule = Constraint(model.S, rule=source_balance_rule)

def flow_balance_rule(model, d):
    return (sum(model.Y[i,d] for i in model.I if model.Arc[i,d]>=1)
            - sum(model.Y[d,i] for i in model.I if model.Arc[d,i]>=1)) == 0.0
model.FlowBalanceRule = Constraint(model.D, rule=flow_balance_rule)

def extra_routing_rule(model,i,j):
    return model.Y[i,j] - model.ArcDemand[i,j] <= model.Z[i,j]
model.ExtraRoutingRule = Constraint(model.I, model.J, rule=extra_routing_rule)

def y_bound_rule(model,i,j):
    return (0.0, model.Y[i,j], model.Arc[i,j] * 51)
model.YBoundRule = Constraint(model.I, model.J, rule=y_bound_rule)


# Stage-specific cost
def first_stage_profit_rule(model):
    return model.FirstStageProfit == 0.0
model.GetFirstStageProfit = Constraint(rule=first_stage_profit_rule)

def second_stage_profit_rule(model):
    return model.SecondStageProfit == \
           sum(sum(model.Rev[i,j] * model.Y[i,j] -
                   (model.Rev[i,j] + model.Cost[i,j])* model.Z[i,j]
                   for i in model.I) for j in model.J)
model.GetSecondStageProfit = Constraint(rule=second_stage_profit_rule)

# Objective
def total_profit_rule(model):
    return (model.FirstStageProfit + model.SecondStageProfit)
model.TotalProfit = Objective(rule=total_profit_rule, sense=maximize)
"""

# Solve WS for given number of sample realizations
numSamples=10
numX = 5

optVal  = numpy.array ([0 for i in range(numSamples)])
optSoln = numpy.array([[0 for i in range(numSamples)] for j in range(numX)])

opt = SolverFactory('gurobi')

for i in range(numSamples):
    datafile = '../scenariodata/Scenario' + str(i+1) + '.dat'
    instance = model.create(datafile)
    results = opt.solve(instance)
    print "Solve" + str(i) + "th instance"
    instance.load(results)
    optVal[i] = value(instance.TotalProfit)
##    for j in range(numX):
#        optSoln[j][i] = value(instance.X[j+1])

# Calculate point est / interval est of objective value, point estimate of solution
z_val = 1.96

WS     = optVal[:].mean()
WS_var = optVal[:].var() * numSamples/(numSamples-1)
WS_halfwidth = z_val*sqrt(WS_var/numSamples)

WS_CI_lo = WS - WS_halfwidth
WS_CI_hi = WS + WS_halfwidth

#X_WS = [0 for j in range(numX)]

#for j in range(numX):
#    X_WS[j] = optSoln[j][:].mean()

print WS
print WS_CI_lo, WS_CI_hi
#print X_WS

